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We have developed an efficient scalable kernel method for thermal transport in open systems, with 
which we have computed the thermal conductance of a junction between bulk silicon and silicon 
nanowires with diameter up to 10 nm. We have devised scaling laws for transmission and reflection 
spectra, which allow us to predict the thermal resistance of bulk-nanowire interfaces with larger cross 
sections than those achievable with atomistic simulations. Our results indicate the characteristic 
size beyond which atomistic systems can be treated accurately by mesoscopic theories. 
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I. INTRODUCTION 

Nanostructures and nanostructured materials offer the 
possibility to tune heat transport properties over an ex- 
ceptionally wide range. For example in carbon based 
materials it is possible to obtain variations of the ther- 
mal transport coefhcients over three orders of magnitude: 
graphene and suspended carbon nanotubes are possibly 
the most efficient heat conductors''^, vifhereas nanotube 
pellets and graphene nanoribbons with disordered edges 
are predicted to have thermal insulating properties^i^. 
Similarly, nanostructuring may turn silicon and SiGe 
alloys into efficient thermoelectric materials, by signifi- 
cantly reducing the thermal conductivity (k) as in the 
case of nanowires^"— (SiNW), SiGe nanocomposites^, su- 
perlatticesS., and nanoporous silico n'"'^' . 

Further improvement in designing materials and nano- 
devices with controlled thermal transport properties 
stems from a deeper theoretical understanding of phonon 
transport. Following Landauer and Biittiker's worksi^ii^, 
atomistic Green's function (GF) formalism has become 
the reference method to study coherent electronic trans- 
porti^"— . The GF approach has been transferred suc- 
cessfully to compute thermal transport in nanostruc- 
turesii""— , and it is the optimal framework to investi- 
gate elastic phonon scattering from impurities, defects, 
disorder or interfaces, i.e. in all those cases where an- 
harmonic phonon-phonon scattering can be deemed of 
secondary importancc'^^'^^. An atomistic GF method in- 
cluding phonon-phonon scattering has also been devel- 
oped and applied to small model systems^^, however one 
can in general safely assume elastic scattering when a 
finite nanoscale system between two reservoirs is con- 
sidered. This is generally the case for molecular junc- 
tions, contacts and grain boundaries. Special care must 
be taken in testing this assumption when one wants to 
extrapolate finite size calculations to extended materials, 
where long wavelength phonons do not get scattered by 
nanoscale impurities and contribute a significant amount 
to the total thermal conductivity. In spite of significant 
insight achieved in these studies based on the GF formal- 
ism, it remains a formidable task to perform atomistic 
simulations of nanostructures with characteristic sizes of 



several tens of nanometers, as it would be needed to 
bridge the gap between theory and experiment. Because 
of matrix inversion operations, even the recursive imple- 
mentation of the GF method, which permits us to deal 
with systems extending for several micrometers in the 
direction of heat propagation, imposes severe size limi- 
tations in the orthogonal plane. In terms of SiNW, this 
means that one is limited to diameters that do not exceed 
few nanometers. Similar limitations hamper the predic- 
tive power of approaches based on molecular dynamics, 
so far restrained to the study of thin wires2ii2^. 

Here we outline a formalism for phonon transport 
based on the scattering-matrix approach^^, which cir- 
cumvents the matrix inversion problem by substitut- 
ing eigenvalue equations with local kernel search and 
intersections. After deriving a generalized scattering- 
matrix approach for phonon propagation, we illustrate 
the numerically stable and efficiently parallelizable kernel 
method. For example, we apply the scattering formalism 
to compute the contact thermal resistance between bulk 
silicon and SiNWs with diameters up to 14 nm. 



II. SCATTERING MATRIX APPROACH 

The scattering-matrix approach was formulated to 
solve quantum electronic transmission problems^S^, and 
found its natural application for the simulation of scan- 
ning tunneling microscopy images^ and of molecular 
electronic devices via the so called elastic scattering 
quantum chemistry (ESQC) method^^. Here we refor- 
mulate the theory in terms of phonon transport. We 
consider a phonon wave packet, represented by a weight- 
normalized displacement field u, traveling through an 
open system made of semi-infinite reservoirs connected 
by an arbitrary structure (defect). Our goal is to deter- 
mine the thermal energy exchanged between the reser- 
voirs through the defect in stationary non-equilibrium 
conditions, i.e. when the reservoirs are kept at different 
temperatures. 

In the harmonic approximation, the equation of motion 
for the displacement field u{t) is il{t) — 'Du{t), where D is 
the dynamical matrix. The real- valued state u can be de- 
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composed in terms of the complex valued eigenstates v{uj) 
of D. Given the state u{tq) and its eigen-decomposition 
coefficients groi^^), the time propagation of u is: 



u{t)= / L„(c.)z;(a.)e— + 



cc. 



dw 



(1) 



Let P be the projector associated with the degrees of 
freedom of an arbitrary part P of the system. To get the 
energy exchanged between P and the rest of the system, 
one can balance the time derivatives of the work from P 
to the whole system and vice versa, thus obtaining: 



Ep{t) = {u{t)\[P,'D]\u{t)) . 



(2) 



The energy of P in stationary conditions {Ep{oo)) is 
found by integrating [2] to the infinite time limit. Substi- 
tuting u with its eigen-decomposition in[l]in the integral 
leads to: 

Ep{^)^-2Tri J naj\gro{io)\^v{u;)\[P,T>]\v{iu))dLu . 

(3) 

All information concerning the initial state lies in the 
weights grg{uj), which can be taken as the statistical dis- 
tribution of the states \v) when simulating a system at fi- 
nite temperature. In the stationary non-equilibrium case, 
those weights refer to the rate of phonons emitted from 
the reservoirs (i.e. ID phonon gas obeying Bose-Einstein 
statistics): 
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where f{uj,T) is the Bose-Einstein distribution function 
at the reservoir temperature T. In order to evaluate [31 
the eigensolutions \v{uj)) of the open system have to be 
expressed in terms of a convenient basis made of a single 
phonon mode \ipl^^{uj)) coming from a reservoir A into 
the defect, and the set of phonon modes ip°'^*{uj) coming 
out of the defect toward the reservoirs: 



(^)iV'r(^)) 



l^'^^(^)), (5) 



where both defect displacements and reservoir surface 



states at the interfaces are included in 



,def 



(lu)). The 



scattering tensor S{uj) maps the incoming phonons 
\ipl'^{uj)) onto the outgoing phonons (uj)) . As the 
energy carried by any incoming or outgoing phonon with 
frequency uj is quantized as [3] provides the following 
normalization and orthogonality conditions: 
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(6) 



where denotes the projector on reservoir A. Com- 
bining the stationary non-equilibrium weights of 2] with 



[3l and observing the conditions of |6l one obtains the sta- 
tionary energy transfer between reservoirs A and B: 

ieAjeB 

(7) 

Once S(ix') is obtained by computing the eigenstates 
w(a;)), the energy flux between two reservoirs A,B is 
determined using the transmission coefficient 7ab(w) = 
'^ieA'^jeB ■ The corresponding thermal con- 

ductance is given by the Landauer formula as the limit 
of El when Ta-^Tb: 
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III. SCALABLE IMPLEMENTATION 




FIG. 1: Silicon nanowire/bulk interconnect partitioning. The 
central defect is generally further subdivided in order to speed 
up the computation. The bulk and wire reservoirs are or- 
ganized as a pile of periodic slices (S and S', respectively) 
indexed starting from the contact areas. 

The first step to compute the eigenstates in [5] is to 
rewrite the eigenvalue problem as a null-space search 
problem: 



Dv = oj'^v <^ w e ker{D - cli^} 



(9) 



It is then possible to consider a partition V — {Pi} of 
the system , and to solve locally the auxiliary equations 



V e ker{Pi(D - w^)}. 



(10) 



The solutions of the eigenvalue problem are given by the 
intersection of the resulting invariant subspaces: 



^ 1- e P|ker{Pi(D -a;2) . (n) 



Typically, the partition is a set of projectors on each 
reservoir, completed by a set of projectors on the defect. 
The auxiliary equations ([TU]) for the defect are solved 
with a QR decomposition method, which is stable and 
numerically efficient. 

The reservoirs are treated in a separate way: as in the 
propagator method^^, every reservoir is partitioned into 
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periodic slices Si of dimension n such that only nearest 
neighboring slices interact (see[T]). However, instead of 
formulating a spatial propagator, we first compute the 
2n non-trivial solutions of [TU] for the second slice 5*2 of 
the reservoir: 



'S2 



-"SI 
'JS2 



^S3 



(12) 



The periodic solutions are then reconstructed by solving 
the 2n x 2n generalized eigen-problem: 
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(13) 

As every slice of the reservoir except SI is equivalent to 
S2, the periodic solutions hold for the entire reservoir, 
except for SI which is treated explicitly as part of the 
defect. The intersection of the periodic solutions leads 
to phonon modes |7/'™/°"*(w)) (|a//3|=l), and surfaces 
states (I a//? I ^ 1). After intersecting the reservoir and 
the defect solutions, we extract a basis {vi\ spanning 
only surface states localized at the defect interface (i.e. 
with |a//3| < 1): 

v. = Y, (A,,|^f (..)) +r,,|^™*(^))) + |?Jf ^(^)). (14) 



The scattering tensor is then computed by applying the 
A~^ transform to the {vi\ set, providing the set of eigen- 
states {vi] defined in[5l so that: S(w) = T ■ A~^. In the 
presence of short range interactions, parts can be defined 
to be as small as the interaction range, so that only neigh- 
boring parts interact. Such an implementation allows 
efhcient parallelization, in the same fashion as domain 
decomposition in molecular dynamics codes. The auxil- 
iary equations are solved locally in parallel before their 
intersection according to a binomial tree. Furthermore, 
a reciprocal space sampling technique allows an efficient 
treatment of the periodic reservoirs solutions. Within 
this framework, the main limitation of the approach is 
the treatment of non-periodic ID reservoirs, which re- 
quires the full diagonalisation of a matrix, which scales 
as the square of the surface of the contact. 



IV. RESULTS AND DISCUSSION 

a. Thermal conductance of bulk silicon/silicon 
nanowire contacts. We apply the scattering matrix 
approach to compute the contact thermal resistance 
of a bulk-SiNW interface. Interface resistance plays 
an essential role in determining the thermal transport 
performance of nanostructured materials and nanoscale 
devices. In addition, evaluating the thermoelectric 
performances of nanostructures such as SiNW, it is 



indispensable to be able to resolve the contact ther- 
mal resistance from the intrinsic resistance. A few 
special cases, such as grain boundaries in silicon, crys- 
talline/amorphous interfaces and silicon/germanium 
junctions, have previously been addressed using 
molecular dynamics and real-space Kubo-Greenwood 
formalis m^^'^° . An often overlooked yet omnipresent 
case where contact resistance is essential is the junction 
between nanostructures and reservoirs. A simplified 
model, based on lattice dynamics calculations of bulk 
silicon and of SiNWs with different diameters, predicts 
that the contact resistance is dominant over the intrinsic 
resistance of the ideal nanowire^i*^. Coherent contacts 
between crystals and nanowires with diameter as small 
as ^ 20 nm can be actually realized by etching nanowires 
directly out of the bulk precurso r^'^i'^^ . We model the 
interatomic interactions between silicon atoms by means 
of the short-range empirical forcefield after Tersof^. 
Crystalline SiNWs with diameters between 2 and 14 
nm are considered. The wires are grown in the (100) 
crystallographic direction, have a nearly circular cross 
section and are coherently connected to the bulk reser- 
voir. The surface is reconstructed in order to minimize 
the number of dangling bonds^S. 

Transmission spectra are displayed in HJa) along with 
the interface conductance a (b) obtained by integrat- 
ing the transmission coefficient over the whole frequency 
spectrum according to HI The data sets are normal- 
ized according to the interface area, as one would rea- 
sonably expect the conductance of a junction to scale 
with its cross section area. In fact such normalization 
makes curves comparable, but not overlapping. Normal- 
ized transmission spectra and ct{T) curves overlap for 
wires of 7 and 10 nm diameter (red and pink curves in[2|). 
Whereas heat transport in thicker wires can be treated 
within a mesoscopic approacb^i^, below this threshold, 
one has to consider explicitly the atomistic details of the 
interface to obtain an accurate estimate of the contact 
conductance. As the construction of the bulk- wire inter- 
face is ideal at the atomic scale, our calculations provide 
an upper limit to the contact conductance. In the low 
temperature regime (T < 50 K) the interface area nor- 
malized contact conductances collapse to a single curve 
and display a temperature dependence of T^. This trend 
was formerly predicted analytically^ and confirmed in 
experiments^^, where it was shown that deviations from 
the behavior stem from specific features of the SiNW, 
such as surface roughness, the effects of which add up in 
series to the contact conductance. On the other hand we 
observe that the reflection spectrum (not shown), scales 
with the linear dimension of the contact interface, which 
seems to indicate that back scattering of phonons mostly 
happens at the perimeter of the junction. Normalized re- 
flection spectra nearly coincide for wires with a diameter 
larger than 7 nm. 

b. Representation of the heat flux. An advantage 
of the present implementation of the scattering matrix 
method is that it provides a real-space representation of 
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FIG. 2: (color) Transmission spectra (a) and thermal con- 
ductance as a function of the temperature (b) for a set of 
bulk-nanowire contacts. Both data sets are normalized with 
respect to the interface area expressed either in nm (conduc- 
tance) or in number of atoms (transmission). 

the energy flux at any given frequency. This allows visu- 
ahzation of the parts of the system that primarily trans- 
mit or reflect thermal energy. An example is shown in 
[3l where the norm of the heat flux across a bulk-lOnm 
SiNW interface is represented. Phonon branches at 0.25, 
0.75, 2 and 4 THz are considered. The spacial features of 
heat transport at different frequencies are clearly differ- 
ent: whereas at the lowest frequency (0.25 THz) thermal 
energy is mainly transmitted through the central bulk- 
like part of the wire, at higher frequencies (0.75 and 2 
THz) thermal energy is transferred through a surface 
layer. Beyond 4 THz heat is transferred through the 
center of the wire. We note that phonons with frequency 
between ~ 1 and 4 THz, which are the majority heat 
carriers in crystalline Si at room temperature, transfer 
energy preferably through a sub-surface layer. Therefore 
our results indicate the reason why thermal conductivity 
of SiNW is so sensitive to surface modifications, such as 
disorder or presence of interfaces^i^. 

c. Dimensionality and shape effects. In order to 
probe the effects of shape and dimensionality reduction 
on the contact conductance we compare the number of 
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FIG. 3: (color) Volumetric representation of the norm of the 
energy flux at the interface of a 10 nm thick silicon nanowire, 
corresponding to channels with frequency of 0.25, 0.75, 2 and 
4 THz. In the 0.75 and 2 THz case, thermal transport mainly 
occurs in a thin sub-surface layer (red color area). 



phonon channels (corresponding to the density of states) 
over the whole frequency spectrum, in contacts made of 
crystalline bulk silicon and either wires with a circular 
section or square rods. We only consider SiNW larger 
than the threshold size of 7 nm, identified as the onset 
for a mesoscopic theory of thermal transport. The cal- 
culations have been performed for SiNW with diameters 
up to 14 nm. The data are conveniently normalized with 
respect to the contact surface area and are compared 
to the number of channels in three-dimensional periodic 
bulk. To verify size convergence we consider two bulk 
samples with cubic supercell of 8.7 and 13 nm, respec- 
tively (|T]). Our data show that for SiNWs larger than 
7 nm, the number of channels per atom at a given fre- 
quency does not depend on the diameter. The number of 
channels at low frequency (< 3 THz) for the contacts is 
the same as in the crystalline bulk, but it deviates signif- 
icantly from the bulk at larger frequencies. This means 
that even in contact interfaces with very large wires, one 
cannot expect to recover bulk-like thermal conductance. 
It also indicates that dimensionality reduction has a pro- 
found effect on the limit density of states as well. Such 
a limit depends also on the shape of the SiNW, but to a 
minor extent. The spectrum of square shaped nanorods 
differs from that of circular ones in the medium-to-high 
frequency range, but it retains similar features as cylin- 
dric wires and does not seem to approach the 3D bulk 
limit either. 



V. CONCLUSIONS 

We have developed an efficient method based on the 
scattering matrix approach to compute the thermal con- 
ductance in an open system. Our derivation leads to 
an expression of the energy flux between two reservoirs 
across a defect region, equivalent to the one derived in 
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traditional Green's function approach, as it circumvents 
the bottleneck imposed by the inversion of large matrices, 
and allows real size devices to be simulated at the atom- 
istic level. We have used this approach to compute the 
contact thermal conductance of ideal junctions between 
bulk silicon and silicon nanowires of different diameters. 
Our results show that beyond a threshold diameter of 7 
nni phonon transmission, reflection and thermal conduc- 
tance obey simple scaling laws, whereas deviations are 
observed for thinner wires. Our approach also provides 
a direct space visualization of frequency dependent heat 
flux, which yields valuable insight into the spatial fea- 
tures of heat conduction in nanoscale devices. 



FIG. 4: (color) Number of transmission channels for a set of 
bulk-nanowire contacts of different diameter and shape. The 
data are normalized with respect to the interface area ex- 
pressed in number of atoms. Data are compared to the num- 
ber of channels in a three-dimensional periodic bulk to high- 
light the effect of dimensionality reduction. All the nanowires 
considered here are larger than the 7 nm diameter threshold. 

RefsJi^— . However, our implementation differs from the 
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